clear all
set more off, permanently


********************************************************************************
**************************** main estimations **********************************
********************************************************************************

*********************************** Roraima ************************************

******* generate dataset that will be used to save Roraima results

clear
set obs 80
egen year=seq(), f(2000) t(2015)

gen n=_n
gen type="."
replace type="C" if n<17
replace type="R" if n>16 & n<33
replace type="U" if n>32 & n<49
replace type="XY_CtoR" if n>48 & n<65
replace type="XY_RtoC" if n>64 & n<=80
drop n

gen C=.
gen R=.
gen U=.
gen XY_CtoR=.
gen XY_RtoC=.

save  "...\roraima_GFC_trend", replace

******* forest cover losses in reduced and constant PAs

clear

use  "...\DID_BZ20km_GFC.dta" 

keep if  ID_state==23 

egen enacR_09_= max(enacR_09), by(ID1)

drop if year<2000

gen PA_=1 if PA!=.
bysort ID1: egen PA0108=max(PA_) if year<2009
bysort ID1: egen PA0108_=max(PA0108) 
drop PA0108
rename PA0108_ PA0108

keep if  PA0108==.

collapse(mean) treecoverloss, by(enacR_09_ year)
separate treecoverloss, by(enacR_09_)

*** prepare data to merge to roraima_GFC_trend

drop treecoverloss

rename enacR_09_ type
tostring type, replace
replace type="C" if type=="0"
replace type="R" if type=="1"
rename treecoverloss0 C
rename treecoverloss1 R
drop if type=="."

merge 1:1 year type using "...\roraima_GFC_trend", nogen

save  "...\roraima_GFC_trend", replace

******* forest cover losses in unprotected areas

clear

use  "...\DID_BZ20km_GFC.dta" 

keep if  ID_state==23

drop if year<2000

collapse(mean) treecoverloss, by(R09protected2000 year)
separate treecoverloss, by(R09protected2000) 

*** prepare data to merge to para_GFC_trend

drop treecoverloss treecoverloss1

rename R09protected2000 type
tostring type, replace
replace type="U" if type=="0"
rename treecoverloss0 U
drop if type=="." | type=="1"

merge 1:1 year type using "...\roraima_GFC_trend", nogen

save  "...\roraima_GFC_trend", replace

******* forest cover losses in matched PAs

clear all

use "...\DIDmatching_BZ20km_GFC.dta" 

egen enacR_09_= max(enacR_09), by(ID1)

keep if  PA0108==. & ID_state==23 

drop if roraima_weight==.

collapse(mean) treecoverloss, by(enacR_09_ year)
separate treecoverloss, by(enacR_09_) 

*** prepare data to merge to roraima_GFC_trend

drop treecoverloss 

rename enacR_09_  type
tostring type, replace
replace type="XY_CtoR" if type=="0"
replace type="XY_RtoC" if type=="1"
rename treecoverloss0 XY_CtoR
rename treecoverloss1 XY_RtoC
drop if type=="." 

merge 1:1 year type using "...\roraima_GFC_trend", nogen

save  "...\roraima_GFC_trend", replace

***** figure

clear

use "...\roraima_GFC_trend"

graph twoway line  C R U XY_CtoR XY_RtoC  year   ///
		, lcolor(gs12 gs8  black gs12 gs8) lpattern(solid solid solid dash dash) ///
		xlabel(2001(1)2015, angle(45)) xtitle("Year") ylabel(0(0.01)0.04, angle(horizontal)) ytitle("Proportion of deforested pixels") ///
		legend(col(2) symxsize(8)  size(vsmall) label(1 "Constant-Sized PAs") label(2 "PA Size Reductions") label(3 "Unprotected Areas") label(4 "Matched Constant-Sized PAs") label(5 "Matched PA Size Reductions") ) ///
		xline(2009, lpattern(dash) lcolor(black)) ///
		graphregion(fcolor(white)) sort
		
graph export "...\figure6A_roraima.png", replace


********************************** Para ***************************************

******* generate dataset that will be used to save Para results

clear
set obs 80
egen year=seq(), f(2000) t(2015)

gen n=_n
gen type="."
replace type="C" if n<17
replace type="R" if n>16 & n<33
replace type="U" if n>32 & n<49
replace type="XY_CtoR" if n>48 & n<65
replace type="XY_RtoC" if n>64 & n<=80
drop n

gen C=.
gen R=.
gen U=.
gen XY_CtoR=.
gen XY_RtoC=.

save  "...\para_GFC_trend", replace

******* forest cover losses in reduced and constant PAs

clear

use  "...\DID_BZ20km_GFC.dta" 

keep if  ID_state==14 

egen enacR_12_= max(enacR_12), by(ID1)

drop if year<2000

gen PA_=1 if PA!=.
bysort ID1: egen PA0108=max(PA_) if year<2009
bysort ID1: egen PA0108_=max(PA0108) 
drop PA0108
rename PA0108_ PA0108

keep if PA0108==.

collapse(mean) treecoverloss, by(enacR_12_ year)
separate treecoverloss, by(enacR_12_)

*** prepare data to merge to para_GFC_trend

drop treecoverloss

rename enacR_12_ type
tostring type, replace
replace type="C" if type=="0"
replace type="R" if type=="1"
rename treecoverloss0 C
rename treecoverloss1 R
drop if type=="."

merge 1:1 year type using "...\para_GFC_trend", nogen

save  "...\para_GFC_trend", replace

******* forest cover losses in unprotected areas

clear

use  "...\DID_BZ20km_GFC.dta" 

keep if  ID_state==14 

drop if year<2000

keep if  R12protected2000!=.  

collapse(mean) treecoverloss, by(R12protected2000 year)
separate treecoverloss, by(R12protected2000)

*** prepare data to merge to para_GFC_trend

drop treecoverloss treecoverloss1

rename R12protected2000 type
tostring type, replace
replace type="U" if type=="0"
rename treecoverloss0 U
drop if type=="." | type=="1"

merge 1:1 year type using "...\para_GFC_trend", nogen

save  "...\para_GFC_trend", replace

******* forest cover losses in matched PAs

clear all

use "...\DIDmatching_BZ20km_GFC.dta" 

egen enacR_12_= max(enacR_12), by(ID1)

keep if  PA0108==. & ID_state==14  

drop if para_weight==.

collapse(mean) treecoverloss, by(enacR_12_ year)
separate treecoverloss, by(enacR_12_) 

*** prepare data to merge to para_GFC_trend

drop treecoverloss 

rename enacR_12_  type
tostring type, replace
replace type="XY_CtoR" if type=="0"
replace type="XY_RtoC" if type=="1"
rename treecoverloss0 XY_CtoR
rename treecoverloss1 XY_RtoC
drop if type=="." 

merge 1:1 year type using "...\para_GFC_trend", nogen

save  "...\para_GFC_trend", replace

*** figure

clear

use "...\para_GFC_trend"

graph twoway line  C R U XY_CtoR XY_RtoC  year   ///
		, lcolor(gs12 gs8  black gs12 gs8) lpattern(solid solid solid dash dash) ///
		xlabel(2001(1)2015, angle(45)) xtitle("Year") ylabel(0(0.01)0.04, angle(horizontal)) ytitle("Proportion of deforested pixels") ///
		legend(col(2) symxsize(8)  size(vsmall) label(1 "Constant-Sized PAs") label(2 "PA Size Reductions") label(3 "Unprotected Areas") label(4 "Matched Constant-Sized PAs") label(5 "Matched PA Size Reductions") ) ///
		xline(2012, lpattern(dash) lcolor(black)) ///
		graphregion(fcolor(white)) sort
		
graph export "...\figure6B_para.png", replace

********************************* Rondonia ************************************* 

******* generate dataset that will be used to save Rondonia results

clear
set obs 80
egen year=seq(), f(2000) t(2015)

gen n=_n
gen type="."
replace type="C" if n<17
replace type="R" if n>16 & n<33
replace type="U" if n>32 & n<49
replace type="XY_CtoR" if n>48 & n<65
replace type="XY_RtoC" if n>64 & n<=80
drop n

gen C=.
gen R=.
gen U=.
gen XY_CtoR=.
gen XY_RtoC=.

save  "...\rondonia_GFC_trend", replace

******* forest cover losses in reduced (state) and constant PAs

clear
use  "...\DID_BZ20km_GFC.dta" 

drop if year<2000
drop if enacR_0912==.

keep if  ID_state==22
keep if PA0108==.

egen enacR_10_= max(enacR_10), by(ID1)

keep if (enacR_10_==1 & state_==1 | enacR_10_==0 & state_==.)

collapse(mean) treecoverloss, by(enacR_10_ year)
separate treecoverloss, by(enacR_10_) 

*** prepare data to merge to rondonia_GFC_trend * we only take 

drop treecoverloss

rename enacR_10_ type
tostring type, replace
replace type="C" if type=="0"
replace type="R" if type=="1"
rename treecoverloss0 C
rename treecoverloss1 R
drop if type=="."

merge 1:1 year type using "...\rondonia_GFC_trend", nogen

save  "...\rondonia_GFC_trend", replace

******* forest cover losses in unprotected areas

clear

use  "...\DID_BZ20km_GFC.dta" 

keep if  ID_state==22 
drop if year<2000

gen R10protected2000=1 if Rprotected2000==1 & year_PADDD==2010
replace R10protected2000=0 if Rprotected2000==0 

collapse(mean) treecoverloss, by(R10protected2000 year)
separate treecoverloss, by(R10protected2000) 

*** prepare data to merge to rondonia_GFC_trend

drop treecoverloss treecoverloss1

rename R10protected2000 type
tostring type, replace
replace type="U" if type=="0"
rename treecoverloss0 U
drop if type=="." | type=="1"

merge 1:1 year type using "...\rondonia_GFC_trend", nogen

save  "...\rondonia_GFC_trend", replace

******* forest cover losses in matched PAs

clear all

use "...\DIDmatching_BZ20km_GFC.dta"

keep if  PA0108==. & ID_state==22

keep if (enacR_10_==1 & state_==1 | enacR_10_==0 & state_==.)

drop if rondonia_weight==.

collapse(mean) treecoverloss, by(enacR_10_ year)
separate treecoverloss, by(enacR_10_) 

*** prepare data to merge to rondonia_GFC_trend

drop treecoverloss 

rename enacR_10_  type
tostring type, replace
replace type="XY_CtoR" if type=="0"
replace type="XY_RtoC" if type=="1"
rename treecoverloss0 XY_CtoR
rename treecoverloss1 XY_RtoC
drop if type=="." 

merge 1:1 year type using "...\rondonia_GFC_trend", nogen

save  "...\rondonia_GFC_trend", replace

*** figure

clear

use "...\rondonia_GFC_trend"

graph twoway line  C R U XY_CtoR XY_RtoC  year   ///
		, lcolor(gs12 gs8  black gs12 gs8) lpattern(solid solid solid dash dash) ///
		xlabel(2001(1)2015, angle(45)) xtitle("Year") ylabel(0(0.01)0.04, angle(horizontal)) ytitle("Proportion of deforested pixels") ///
		legend(col(2) symxsize(8)  size(vsmall) label(1 "Constant-Sized PAs") label(2 "State PA Size Reductions") label(3 "Unprotected Areas") label(4 "Matched Constant-Sized PAs") label(5 "Matched PA Size Reductions") ) ///
		xline(2010, lpattern(dash) lcolor(black)) ///
		graphregion(fcolor(white)) sort
		
graph export  "...\figure6C_rondonia.png", replace



